function updateIB_AD(timeStep)
%updateIB_BD solves the diffEq for t = timeStep+dt and returns it
%   Detailed explanation goes here

global t;
global LTSCells;
global IB_ADCells;
global IB_SOCells;
global Synapses;
global LTSIndices;
global IBIndices;
global Gates;
global J_IB;

dt = t(timeStep+1) - t(timeStep);

for i = 1:size(IB_ADCells,3)
    in_s = Synapses(timeStep, :, IBIndices(i));
    in_gates = Gates(:, IBIndices)';
    V_c = IB_SOCells(timeStep,1,i);
    [k1, s1] = feval(@IB_AD, timeStep, t(timeStep), IB_ADCells(timeStep,:,i), in_s, in_gates, V_c, J_IB);
    % note - caromk 5/24/2010 - these k and s updates are to normalize the
    % RK4 method with the cleaned up code, this code should not continue
    % to be built upon
    k1 = dt*k1; s1 = dt*s1;
    [k2, s2] = feval(@IB_AD, timeStep, t(timeStep) + dt/2, IB_ADCells(timeStep,:,i) + (k1)'/2,in_s+s1/2, in_gates, V_c, J_IB);
    k2 = dt*k2; s2 = dt*s2;
    [k3, s3] = feval(@IB_AD, timeStep, t(timeStep) + dt/2, IB_ADCells(timeStep,:,i)+ (k2)'/2,in_s+s2/2, in_gates, V_c, J_IB);
    k3 = dt*k3; s3 = dt*s3;
    [k4, s4] = feval(@IB_AD, timeStep, t(timeStep) + dt, IB_ADCells(timeStep,:,i) + k3',in_s+s3, in_gates, V_c, J_IB);
    k4 = dt*k4; s4 = dt*s4;
    IB_ADCells(timeStep+1,:,i) = IB_ADCells(timeStep,:,i) + (k1' + 2*k2' + 2*k3' + k4')/6;
    Synapses(timeStep+1,:,IBIndices(i)) = Synapses(timeStep,:,IBIndices(i)) + (s1 + 2*s2 + 2*s3 + s4)/6;
end

end